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Abstract —  In  this  paper,  the  problem  of  estimating  the 
frequencies,  dampings,  amplitudes  and  phases  of  closely 
spaced  complex  damped  exponentials  in  the  presence  of 
noise  is  considered.  In  several  papers,  decimation  is  pro¬ 
posed  as  a  way  to  increase  the  performance  of  subspace- 
based  parameter  estimation  methods,  in  the  case  of  over- 
sampling  [1]  [2]  [3] .  In  this  paper,  a  novel  extension  of  the 
HTLS-method  [4]  that  operates  directly  on  the  decimated 
data  matrix  is  presented,  and  it  is  compared  to  other  dec¬ 
imation  methods.  Experiments  on  simulated  nuclear  mag¬ 
netic  resonance  (NMR)  spectroscopy  signals  show  the  in¬ 
fluence  of  decimation  on  the  accuracy  and  computational 
complexity  of  the  estimators. 

Keywords —  decimation,  subspace-based  parameter  esti¬ 
mation,  NMR,  quantitation 


I.  Introduction 

In  various  applications  of  digital  signal  processing,  such 
as  NMR  spectroscopy  and  speech  processing,  complex 
damped  exponentials  are  used  as  a  model  function.  Let 
x(t)  be  a  sum  of  K  complex  damped  exponentials,  con¬ 
taminated  by  additive  white  noise  n(t): 

k 

*(t)  =  5>4  +  »W,  t  =  0,1,...,  N-l  (1) 

k= 1 

with  complex  amplitudes  c^,  k  =  1, . . . ,  K 

Ck  =  afeeJV°+^)  (2) 

and  signal  poles  z/,,  k  =  1, . . .  ,K 

Zk  =  e«2irA-<fc)//.ampl«  (3) 

where  represents  the  amplitude,  (<po  +  <j>k)  the  phase, 
fk  the  frequency  and  dk  the  damping  of  the  kt/l  compo¬ 
nent,  and  / sample  is  the  sampling  frequency,  (fio  is  the 
zero  order  phase,  whereas  (j>k  represents  extra  degrees  of 
freedom  that  may  be  required  under  certain  experimen¬ 
tal  conditions  (usually  all  </>&  are  zero).  The  problem  is 
to  estimate  these  parameters  given  a  set  of  N  noisy  data 
points  x(t),  t  =  0, 1, . . . ,  N—  1. 

It  is  known  that  subspace-based  parameter  estimation 
techniques  perform  poorly  when  applied  to  a  signal  con¬ 
sisting  of  a  sum  of  closely  spaced  complex  damped  ex¬ 
ponentials  [5].  Therefore,  in  recent  publications  different 
decimative  approaches  were  proposed  in  order  to  increase 
the  performance  of  these  subspace-based  methods.  The 
idea  is  to  artificially  increase  the  frequency  separation  by 
decimating  (downsampling)  the  signal,  however  making 
sure  that  no  aliasing  is  introduced. 

In  this  paper  a  novel  extension  of  the  HTLS-method  [4] 
that  operates  directly  on  the  decimated  data-matrix  is 
presented,  and  it  is  compared  to  existing  decimative 
subspace-based  algorithms  [1]  [2]  [3] .  Furthermore,  the  in¬ 
fluence  of  decimation  on  the  accuracy  and  computational 


complexity  of  these  subspace-based  estimators  is  ana¬ 
lyzed.  Extensive  Monte-Carlo  simulations  on  simulated 
NMR  signals  show  the  benefits  of  decimation  in  the  field 
of  NMR  spectroscopy. 


II.  Decimative  methods 
Three  subspace-based  methods  are  briefly  described: 

•  ETLSD:  the  ESPRIT-Total  Least  Squares  algorithm  [6] 
applied  to  decimated  data,  as  presented  in  [1]  [2]; 

•  HTLSD:  the  novel  extension  of  the  Hankel- Total  Least 
Squares  algorithm  [4]  for  decimated  data; 

•  DESE:  another  decimative  subspace-based  parameter 
estimation  algorithm,  recently  proposed  as  Decimative 
Spectral  Estimation  [3]. 

In  what  follows,  scalars  are  represented  by  lower-case 
letters,  vectors  by  bold  lower  case  letters  and  matrices  by 
bold  uppercase  letters.  Furthermore,  a  Matlab  like  nota¬ 
tion  is  used:  a(i)  stands  for  the  \th  element  of  vector  a. 


A.  ETLSD 

The  approach  described  in  [1]  uses  several  decimated  se¬ 
quences  to  calculate  the  sample  covariance  matrix,  which 
is  used  to  estimate  the  frequencies  and  dampings  by 
means  of  the  ESPRIT-TLS  method  [6] . 

The  original  data  sequence  x=  [x(0),  x(l) . . . ,  x(N— 1)], 
can  be  divided  into  D  different  decimated  sequences 
Xi  G  *  =  0, 1, . . . ,  D  —  l: 

x;  =  [x(i),x(D  +  *),...,  x((N/D  —  1  )D  +  *)]  (4) 

with  D  the  decimation  factor,  which  should  be  chosen 
such  that  \fk\  <  (/ sample) /(2Z>),  j  =  1, 2, . . . ,  K  in  order 
to  avoid  aliasing.  Using  the  model  (1),  we  can  write: 

K 

*,(<)  =  £  Ckz[Dt+l>)  +n(Dt  +  i),  t  =  0, 1, . . .  ,N/D  —  1 

k= l 
K 

=  ^2(cki){z'ky  +  n(Dt  +  i)  (5) 

k= 1 


where  Cki  =  CkZlk,  k  =  1, 2, . . . ,  K,  i  =  0, 1, . . . ,  D—  1  and 
z'k  =  zk  ,  k  =  1, 2, . . . ,  K.  From  each  of  these  decimated 
sequences  a  sample  covariance  matrix  with  m  lags  is 
formed,  which  can  be  written  as  the  product  of  Hankel- 
matrices:  R;  =  X;Xf  (6) 

with 


2k(0)  Xi(  1) 

Xi( 1)  X{( 2) 


Xi(N/D—m) 
Xi(N/D—m+ 1) 


Xi(m—  1)  Xi(m) 


Xi(N/D-l)  \ 

(7) 
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where  X?  represents  the  conjugate  transpose  of  Xi.  The 
averaged  sample  covariance  matrix 

1  0-1 

R=nER‘  I8) 

i=0 

is  then  used  to  estimate  the  decimated  signal  poles 
z'k.  k  =  1,  2, . . . ,  K  by  means  of  the  ESPRIT-TLS  algo¬ 
rithm  [6].  From  these  estimates  z'k,  the  estimates  of  the 
original  poles  zk,  and  hence  estimates  of  the  frequencies 
and  dampings,  are  easily  obtained.  The  phases  and  the 
amplitudes,  contained  in  ck,  are  then  calculated  as  the 
least  squares  solution  to  (1),  with  zk  replaced  by  the  es¬ 
timates  Zk- 

Since  the  main  computational  cost  of  ETLSD  is  the 
eigenvalue  decomposition  of  the  m  x  m  sample  covariance 
matrix  R,  it  is  clear  that  the  computational  complexity 
mainly  depends  on  the  parameter  m. 

The  method  presented  in  [2]  is  identical  using  the  same 
covariance  matrix,  calculated  in  a  slightly  different  way. 


B.  HTLSD 

A  similar  subspace-based  method  that  operates  directly 
on  the  data  matrix  is  HTLS  [4].  The  difference  between 
HTLS  and  ESPRIT-TLS  is  that  HTLS  makes  use  of  the 
singular  value  decomposition  (SVD)  of  the  data  matrix 
instead  of  the  eigendecomposition  of  the  sample  covari¬ 
ance  matrix. 

For  noiseless  data,  using  (5) ,  X,  can  be  written  in  terms 
of  Vandermonde  matrices: 


Xi 


1  1  ...  1 

A  A  •  •  •  z'k 


Z 


/  m—  1 
1 


Z 


/  m—  1 
2 


CK 


m—  1 


Cn  0 

0  c2i 

0  0 


SC,Tt 


- 

'l 

A  • 

,  N/D-m  “ 
•  Z1 

1 

A  ■ 

/  N/D—m 
•  z2 

1 

z'k  ■ 

,  N/D-m 
■  ZK 

(9) 


block-Hankel  matrix  ~Kstack,  constructed  as  follows: 

xstacfc  =  [  Xi  X2  ...  XD  ]  (10) 

=  S[CXT  C2T  ...  CdT],  (11) 


Since  in  this  case  all  samples  are  used  for  the  SVD,  the 
estimated  poles  are  more  accurate  than  those  estimated 
from  only  one  decimated  sequence. 

After  the  estimation  of  the  signal  poles,  the  phases  and 
the  amplitudes  are  calculated  as  the  least  squares  solution 
to  (1),  with  Zk  replaced  by  the  estimates  Zk- 

Since  the  main  computational  cost  of  HTLSD  is  the 
SVD  of  the  mx  ( N  +  D  —  mD )  data  matrix  X stack,  it 
is  clear  that  the  computational  complexity,  for  given  N, 
mainly  depends  on  the  parameters  m  and  D. 

The  fact  that  HTLS  applied  to  Xstacfc  is  the  decimative 
version  of  HTLS  corresponding  to  ETLSD,  can  also  be 
deduced  from  the  observation  that  the  averaged  sample 
covariance  matrix  R  (8)  can  be  written  as  follows: 


R=^[Xi  X2  ...  XD] 
C.  DESE 


X^ 

Xi 


—  jj^-stack^- 


* 

stack ' 


X 


D  J 


(12) 


This  algorithm  was  presented  very  recently  [3].  Like 
HTLS,  DESE  also  makes  use  of  the  SVD  of  a  Hankel  ma¬ 
trix  and  the  full  set  of  data. 

A  Hankel  matrix  X  is  constructed  from  the  original 
data  sequence  (as  in  (7)  with  D  —  1).  From  this  Hankel 
matrix,  XD  and  Xd  are  computed  by  deleting  respec¬ 
tively  the  top  and  bottom  D  rows.  DESE  uses  the  shift- 
invariance  between  XD  and  Xd  in  order  to  estimate  the 
decimated  poles. 

Without  decimation  (D  =  1) ,  this  method  is  identical 
to  a  method  called  MATPEN,  proposed  in  [7]. 

The  main  computational  cost  of  DESE  consists  of  the 
SVD  of  the  (m—D)x(N—m+ 1)  data  matrix  Xd  and  the 
eigendecomposition  of  an  ( m—D )  x  ( m—D )  matrix.  Since 
D  is  usually  significantly  smaller  than  m,  DESE  usually 
requires  more  operations  than  ETLSD  or  HTLSD. 


From  this  Vandermonde  decomposition,  the  decimated 
signal  poles  z'k,  k  =  1, K  and  complex  amplitudes 
Cki ,  k  =  1  can  immediately  be  derived.  How¬ 

ever,  no  algorithm  exists  to  compute  the  Vandermonde 
decomposition  directly.  Therefore,  HTLS  makes  use  of 
the  shift-invariant  structure  of  S  and  the  singular  value 
decomposition  of  Xi  to  determine  the  decimated  pole  es¬ 
timates  z'k,  k  =  1, 2, . . . ,  K  [4], 

From  (9)  it  is  clear  that  the  Vandermonde  decompo¬ 
sition  of  every  Xi,  *  =  1,2,  has  the  same  (shift- 

invariant)  matrix  S,  i.e.  all  decimated  sequences  have  the 
same  poles.  Therefore,  HTLS  can  also  be  applied  to  the 


III.  Experimental  results 

Extensive  simulations  have  been  performed  on  typical 
simulated  NMR  signals.  Below,  one  representative  exam¬ 
ple  simulating  a  typical  5  peak  31P  NMR  signal  of  per¬ 
fused  rat  liver,  is  given.  N  data  points  (here,  AT  =  128), 
uniformly  sampled  at  10  kHz,  are  generated  by  a  fifth  or¬ 
der  (K=5)  model  function  (1),  of  which  the  parameters 
are  displayed  in  Table  1.  For  several  combinations  of  noise 
level  er„,  matrix  dimension  in  and  decimation  factor  D , 
the  three  described  methods  are  compared  by  means  of 
Monte-Carlo  simulations  consisting  of  2000  noise  realiza¬ 
tions  each.  The  data  points  are  perturbed  by  Gaussian 
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TABLE  I 

Exact  parameter  values  of  the  simulated  NMR  signal 


peak 

/fc(Hz) 

dk  (rad/s) 

ofc(a.u.)a 

M°r 

1 

-1379 

208 

6.1 

15 

2 

-685 

256 

9.9 

15 

3 

-271 

197 

6.0 

15 

4 

353 

117 

2.8 

15 

5 

478 

808 

17.0 

15 

aa.u.  means  arbitrary  units. 

b/ipk  =  </>o  *  180 / 7r  expresses  the  phase  in  degrees 
(in  this  example  </>&  =  0,  k  =  1,  .  .  .  ,  K ) . 

noise  whose  real  and  imaginary  components  have  stan¬ 
dard  deviation  oy.  (Relative)  root  mean  squared  errors 
of  the  estimates  of  all  signal  parameters  are  calculated 
as  well  as  the  percentage  of  failures  per  noise  level.  A 
failure  occurs  when  not  all  5  peaks  are  resolved  within 
specified  intervals  lying  symmetrically  around  the  exact 
frequencies,  or  when  the  estimated  damping  is  negative. 
The  halfwidths  of  the  intervals  are  based  on  Cramer- Rao 
lower  bound  considerations,  and  are  respectively  82,  82, 
82,  43  and  82  Hz. 

The  root  mean  squared  error  (RMSE)  of  the  frequency 
estimates  (excluding  failures)  of  peak  4  is  plotted  as  a 
function  of  ay  in  Fig.  1,  for  ETLSD  and  HTLSD  and 
for  different  values  of  in  and  D.  From  the  curves  with 
m  =  32  it  can  be  seen  that  decimation  increases  the  sta¬ 
tistical  performance  significantly,  both  for  ETLSD  and 
HTLSD.  There  is  no  significant  difference  in  accuracy  be¬ 
tween  ETLSD  and  HTLSD.  The  advantage  of  HTLSD 
over  ETLSD  is  that  squaring  the  data,  and  the  associ¬ 
ated  numerical  problems  with  ill-conditioned  matrices,  is 
avoided.  Keeping  D  constant,  the  accuracy  of  the  fre¬ 
quency  estimates  increases  with  increasing  m,  and  hence 
increasing  computational  burden  (compare  the  curves 
with  to  =  32,  D  =  1  to  those  with  to  =  64,  D  =  1). 


Fig.  1.  Plot  of  RMSE  of  the  frequency  estimates  of  peak  4  versus  noise 
standard  deviation  <ju,  for  ETLSD  (A)  and  HTLSD  (x),  for  different 
decimation  factors  D  and  matrix  dimensions  m  (D  =  l,m  =  32:  dash- 
dotted  line;  D  =  l,m  =  64:  solid  line;  D  =  2,  m  =  32:  dashed  line). 


Comparing  the  curves  with  to  =  64,  D  =  1  to  those  with 
to  =  32,  D  —  2,  we  see  that  with  a  lower  to  a  similar  statis¬ 
tical  accuracy  can  be  obtained  by  proportionally  increas¬ 
ing  D.  Simulations  with  varying  to  and  fixed  D  show 
that,  for  ETLSD  and  HTLSD,  the  most  accurate  esti¬ 
mates  are  obtained  for  m  =  i.e.  for  square  data  ma¬ 
trices  Xy  In  this  case,  the  statistical  accuracy  of  the  esti¬ 
mates  obtained  with  or  without  decimation  is  comparable. 
However,  decimation  allows  to  obtain  these  estimates  at 
a  much  lower  computational  cost.  Indeed,  for  to  = 
the  computational  complexity  is  mainly  determined  by 
the  eigendecomposition  of  the  (jjj)  x  (jjj)  matrix  R  for 
ETLSD,  and  by  the  SVD  of  the  (jjj)  x  (y  +  £))  matrix 
Xstack  for  HTLSD.  So,  for  increasing  D,  the  computa¬ 
tional  cost  to  obtain  the  best  possible  estimates  decreases. 
Therefore,  choosing  maximal  D  (in  order  to  avoid  aliasing 
D  <  (f sample) /(2|/fc|),  k  =  1, . . . ,  K)  gives  optimal  perfor¬ 
mance,  for  ETLSD  and  HTLSD,  both  in  statistical  and 
computational  sense. 

With  to=  the  parameter  estimates  obtained  with 
DESE  are  much  less  accurate,  as  can  be  seen  in  Fig.  2. 
In  order  to  obtain  a  comparable  accuracy  with  DESE, 
the  matrices  XD  and  Xd  should  be  as  square  as  possible, 
i.e.  to  =  N^D  ■  However,  since  usually  D«N,  increas¬ 
ing  D  hardly  reduces  the  computational  cost  because  this 
cost  mainly  consists  of  the  calculation  of  the  SVD  of  the 
( iv~D)  x  ( N~D )  matrix  XD  and  the  eigendecomposition 
of  an  (to  —  D)  x  (to  —  D)  matrix.  Using  m  =  ^  on  the 
other  hand  (as  for  ETLSD  and  HTLSD),  lowers  the  com¬ 
putational  burden  but  then  the  accuracy  of  the  parameter 
estimates  decreases  drastically,  as  shown  in  Fig.  2. 

Some  simple  tests  using  the  counter  of  floating  points 
operations  (flops)  in  Matlab  illustrate  the  dependence  of 
the  computational  complexity  of  each  of  the  algorithms 
on  D  and  to.  ,  as  shown  in  Table  II.  Although  more  effi- 


Fig.  2.  Plot  of  RMSE  of  the  frequency  estimates  of  peak  4  versus 
noise  standard  deviation  <ru,  for  HTLSD  (x)  and  DESE  (o),  for  different 
decimation  factors  D  and  matrix  dimensions  m  ( D  —  l,m  =  32:  dash- 
dotted  line;  D  =  l,m  =  64:  solid  line;  D  =  2,m  =  32:  dashed  line; 
D  =  2,m  =  64:  dotted  line  (for  DESE  only)). 
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TABLE  II 

Computational  complexity  of  the  algorithms,  for 

DIFFERENT  D  AND  m  VALUES,  EXPRESSED  AS  FLOPS/106. 


D 

m 

ETLSD 

HTLSD 

DESE 

1 

32 

3.9 

6.8 

13.9 

2 

32 

3.6 

3.9 

13.9 

1 

64 

24.2 

12.2 

91.1 

2 

64 

n.a.a 

n.a.a 

91.0 

an.a.  means  not  available:  for  ETLSD  and 
HTLSD,  m  must  be  smaller  than  N/D. 


cient  implementations  are  possible  for  each  algorithm  (ex¬ 
ploiting  the  Hankel  matrix  structure,  using  partial  SVD 
algorithms),  the  effect  of  varying  D  and/or  m  on  the  com¬ 
putational  burden  will  qualitatively  be  the  same  as  indi¬ 
cated  in  Table  II. 

On  the  other  hand,  DESE  has  less  failures  than  HTLSD 
(and  ETLSD,  which  have  a  comparable  percentage  of  fail¬ 
ures),  as  shown  in  Fig.  3.  For  DESE,  the  number  of 
failures  increases  with  decreasing  D;  for  ETLSD/HTLSD 
the  number  of  failures  decreases  slightly  with  decreasing  D 
(the  differences  are  small).  However,  DESE  without  dec¬ 
imation  is  still  more  robust  than  ETLSD/HTLSD  with 
decimation  (for  constant  m). 

IV.  Conclusions 

In  this  paper,  we  presented  a  decimative  extension  to 
the  HTLS  method  [4],  and  compared  it  with  two  other, 
recently  proposed,  decimative  subspace-based  parameter 
estimation  methods:  ETLSD  [1]  [2]  and  DESE  [3].  The 
principles  of  the  different  methods  are  presented  within 
the  same  framework,  and  the  relations  between  the  meth¬ 
ods  are  pointed  out.  The  algorithms  were  tested  by  means 
of  various  Monte-Carlo  simulations.  It  is  demonstrated 
that  the  dimensions  of  the  data  matrix  (the  number  of 


Fig.  3.  Plot  of  percentage  of  failures  versus  noise  standard  deviation 
cTjy ,  for  HTLSD  (x)  and  DESE  (o),  for  different  decimation  factors  D  and 
matrix  dimensions  m  ( D  —  1,  m  =  32:  dash-dotted  line;  D  =  2,  m  =  32: 
dashed  line;  D  =  1,  m  =  64:  solid  line;  D  =  2,  m  =  64:  dotted  line  (for 
DESE  only)). 


lags  in  the  covariance  matrix)  is  by  far  the  most  important 
factor  determining  the  statistical  accuracy  of  the  param¬ 
eter  estimates,  leading  to  the  conclusion  that  square  data 
matrices  give  the  best  parameter  estimates.  From  this 
observation,  it  can  be  derived  that  decimation  as  applied 
in  ETLSD  and  HTLSD,  does  not  lead  to  better  estimates 
than  the  best  possible  ones  without  decimation.  However, 
decimation  allows  to  obtain  estimates  with  the  same  accu¬ 
racy  at  much  lower  computational  cost.  The  decimative 
approach  DESE,  presented  in  [3],  however,  does  not  have 
this  computational  advantage,  and  its  statistical  accuracy 
did  not  prove  to  be  higher  than  that  of  the  other  two  ap¬ 
proaches.  For  very  high  noise  levels  the  number  of  fail¬ 
ures  for  DESE  seems  to  be  lower  than  that  of  ETLSD  and 
HTLSD.  It  is  however  not  clear  how  any  of  these  methods 
can  be  used  in  practice  (e.g.  NMR  quantitation)  when  the 
percentage  of  failures  of  the  method  is  larger  than  10%, 
meaning  that  one  out  of  ten  quantitations  is  useless. 
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